# Authors:  Neil Malhotra and Jake Jares
# Date:     2021.09.10
# Inputs:   "Survey Respondent Data for APSR Analyses.dta"
# Outputs: "Figure 2.png"
rm(list=ls())
library(foreign)
library(haven)
library(tidyverse)
library(estimatr)

# CHANGE THE FOLLOWING FILEPATH TO REFLECT LOCATION OF REPLICATION FOLDER
PROJECT_DIR = "/Users/jjares/Documents/research/Farm Subsidies/replication_materials"
regression_data_fpath = paste0(PROJECT_DIR, "/data/Survey Respondent Data for APSR Analyses.dta")

# Read in regression dataset, and drop observations that would be dropped
# by the regression anyways
keep_fields = c(
    "complete",
    "mfp_support", "mfp_any", "conservative", "conservative_mfp_any", "military",
    "gender", "age", "education", "total_acres", "farm_value_2019"
)
df = read_dta(regression_data_fpath) %>%
     dplyr::select(all_of(keep_fields)) %>%
     filter(complete == 1) %>%
     drop_na()

# Regress MFP support ~ (binary) MFP reciept
# Adjust standard errors for heteroskedasticity (HC1, same as Stata robust standard errors)
mod = lm_robust(mfp_support ~ mfp_any + conservative + conservative_mfp_any + military + 
                       gender + age + education + total_acres + farm_value_2019,
         data=df, se_type="HC1")

subgroups_with_means = df %>%
    group_by(mfp_any, conservative) %>%
    summarise_all(.funs = c(mean="mean")) %>%
    rename_with(~str_remove(., '.mean'))
    
# Predict outcomes for (MFP, No MFP) X (Conservative, Not Conservative),
# with all other fields set at population means
# NOTE: this is equivalent to margins, at(mfp_any=0 conservative=0) in Stata
subgroups_with_pop_means = subgroups_with_means %>%
                           dplyr::select(mfp_any, conservative, conservative_mfp_any, mfp_support) %>%
                           mutate(military=mean(df$military),
                                  gender=mean(df$gender),
                                  age=mean(df$age),
                                  education=mean(df$education),
                                  total_acres=mean(df$total_acres),
                                  farm_value_2019=mean(df$farm_value_2019))

mfp_predictions = predict(mod, newdata=subgroups_with_pop_means, interval="confidence") %>%
                  data.frame() %>%
                  rename(mfp_support=fit.fit) %>%
                  mutate(mfp_any = subgroups_with_pop_means$mfp_any,
                         conservative = subgroups_with_pop_means$conservative,
                         mfp = factor(mfp_any, levels=c(0,1), labels=c("Did Not Receive MFP", "Received MFP")),
                         ideology = factor(conservative, levels=c(0,1), labels=c("Moderates/Liberals","Conservatives")))

f2 = ggplot(data=mfp_predictions, 
            aes(x=mfp, y=mfp_support, ymin=fit.lwr, ymax=fit.upr, fill=ideology,
                shape=ideology, color=ideology, width = 0.5)) +
    geom_pointrange(size=0.5, position = position_dodge(width=0.2)) +
    theme(text = element_text(size = 16)) +
    scale_color_manual(values=c("black","black")) +
    scale_x_discrete(name="") +
    scale_shape_manual(values=c(21,24)) +
    scale_fill_manual(values=c("black","white")) +
    ylab("Support for MFP") +
    labs(fill="Ideology",shape="Ideology",color="Ideology")

ggsave(paste0(PROJECT_DIR, "/graphics/Figure 2.tiff"),
       dpi=1000,
       plot=f2)
ggsave(paste0(PROJECT_DIR, "/graphics/Figure 2.pdf"),
       plot=f2)



